Mon. Not. R. Astron. Soc. OOP. [TlfTil (20081 Printed 18 May 2009 



(MN IMtjX style file v2.2) 



Dynamic Monte Carlo radiation transfer in SPH. 
Radiation pressure force implementation. 

Sergei Nayakshin*, Seung-Hoon Cha, and Alexander Hobbs 

Department of Physics & Astronomy, University of Leicester, Leicester, LEI 7RH, UK 



Received 



ABSTRACT 

We present a new framework for radiation hydrodynamics simulations. Gas dynam- 
ics is modelled by the Smoothed Particle Hydrodynamics (SPH), whereas radiation 
transfer is simulated via a time-dependent Monte-Carlo approach that traces photon 
packets. As a first step in the development of the method, in this paper we consider the 
momentum transfer between radiation field and gas, which is important for systems 
where radiation pressure is high. There is no fundamental limitations on the num- 
ber of radiation sources, geometry or the optical depth of the problems that can be 
studied with the method. However, as expected for any Monte-Carlo transfer scheme, 
stochastic noise presents a serious limitation. We present a number of tests that show 
that the errors of the method can be estimated accurately by considering Poisson 
noise fluctuations in the number of photon packets that SPH particles interact with 
per dynamical time. It is found that for a reasonable accuracy the momentum carried 
by photon packets must be much smaller than a typical momentum of SPH particles. 
We discuss numerical limitations of the code, and future steps that can be taken to 
improve performance and applicability of the method. 

Key words: Physical Data and Processes: radiative transfer - Physical Data and 
Processes: hydrodynamics 



1 INTRODUCTION 

Astronomers do not have a luxury of testing their ideas of, 
e.g., galaxy and star formation in a purpose designed lab- 
oratory. Instead, testing grounds are provided by observa- 
tions and numerical simulations. The latter in effect repre- 
sent experiments with a given set of physical laws included. 
Clearly, the more physics is included in the simulations the 
more realistic the latter should be. Presently, two basic phys- 
ical processes - gravity and hydrodynamics - are modelled 
well by a variety of methods and codes. For example, one 
reliable and widely used numerical method of modelling 
gas dynamics is Smoothed Particle Hyd rodynamics (SPH) 
iGingold fc Monaghan|[l977l : ILucvII1977| ). It has been used 
in the various fields of astrophysics, and is especially power- 
ful when resolving high density regions is key. The method 
is grid-less and fully Lagrangian in nature, facilitating mod- 
elling of arbitrary geometr y systems. Fo r reviews of the SPH 
method see , for example , iBenj (|l990t ): iMonaghanI (|l992h : 
iFulkl (|l994 ): |Pricj l|2004 ). 

Radiation transfer and interaction with matter is an- 
other basic process operating in ast rophy si cal s ystems. 
Starting from the pioneering ideas of ILucvI (|l97'i1 ). great 
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efforts have been expended to model radiation-matter 
interactions in SPH. When gaseous systems modelled 
are optically thick, photons scatter or get absorbed and 
rc-cmittcd multiple times before escaping the system 
( R ybicki fc Lig htmaiil Il986l ) . This situation is well approx- 
imated by the diffusion approximation. Several authors 
have alr eady incorporated this radia t ion transfer scheme 
in SPH (lLucvlll977l: lBrookshawlll985l: IWhitehouse fc Bate! 
" I2OO5I : K 
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Viau et al.ll2006h . The diffusion approximation is usually 



complemented by a flux limiter method to model op- 
tically thin regio ns where photons stream freely rather 
than diffuse (e.g., Whitehouse fc Bat jbOO^ ). Most recently 
iPetkova fc Springell ( 2008l 'l implemented the diffusion radia - 
tion scheme in the cosmological code Gadget (|Springelll2005l ) 
supplementing it by the variable Eddington tensor as a clo- 
sure relation. 

Ray tracing methods is a principally different approach 
to radiative transfer, where the radiative transfer equation is 
solved along chosen directions (rays) . Most applications only 
consider rays that start at discrete sources of radiation field, 
such as bright stars, etc., essentially neglecting the diffuse 
radiation field. As far as the density estimation is concerned, 
some authors used density defin ed at the locations of the 
nearest neighbours along the rays (|Kessel-Devnet fc Burkertl 
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I2OOOI : iDale et alj l2007l : iGritschneder et alj 120081 ) . whereas 
others approximated the density field usins the nodes of 
the tree structure one way or another JOxlev fc WoolfsonI 
l2003l :[ Stamatellos fc Whitworthll2005l : ISusall2006l ). The most 
recent developments use the SPH densit y field directly, 
e.g.. c.f. codes SPHRAY dAltav et al.ll2008l) a i id TR APHIC 
(|Pawhk fc Schavd I2OO8I ). iRitzerveld fc IcQ (|2006l ) devel- 
oped a method to transport radiation on adaptive random 
lattices, and showed that the algorith is very efficient for 
cosmological re-ionisation problems. 

Monte-Carlo methods are similar in spirit and yet sub- 
stantially different from the ray tracing codes. In the for- 
mer, the idea is to discretise the radiation field into "pack- 
ets", choose directions and emission time of these packets 
stochastically to obey proper physical constraints, and then 
propagate the packets through matter in accord with radi- 
ation transfer equations. The main problem for the method 
is stochastic noise caused by photon packet statisti cs (for an 
early SPH application see lLucvlll999l ). [BaesI (|2008h presents 
some new interesting ideas about using a smoothing kernel 
in a Monte-Carlo simulation. 

Most of these applications were tailored to photo- 
ionisation problems in the field of star formation or cos- 
mology, where the density field can be considered static as 
rays propagate throu gh it (the "static diffusion limit"; see 
iKrumholz et al.ll2007l '). and where radiation pressure effects 
can be omitted. There are also a number of radiation trans- 
fer codes workin g with a grid-ba sed hydrodynamics codes 
rather SPH (e.g.. llliev et"aLll2006h . We do not discuss these 
methods here. 

Here we present and test a new photon packet based 
radiation transfer scheme combined with an SPH code. The 
method uses the SPH density field directly, thus preserving 
the "native" SPH resolution. The previous efforts did not 
consider the radiation pressure effects, studying instead ra- 
diative heating/cooling and photo-ionisation. Here we study 
the radiation pressure effects instead. We assume that gas 
equation of state is known (e.g., isothermal or polytrophic) 
and consider only the radiation pressure forces. This allows 
us to thoroughly test precision of our approach and numer- 
ical noise effects. There is no fundamental difficulty in in- 
cluding the heating/cooling and photoionisation processes 
in our scheme, and we shall extend our method in that di- 
rection in the near future. Another defining characteristic of 
our new method is that the photon field is evolved in the 
same time-dependent way, although on shorter time steps, 
as gas dynamics, and therefore the method is intrinsically 
time-dependent . 



2 DESCRIPTION OF THE METHOD 
2.1 Radiation Transfer Method 

Radiative transfer equation along a ray 
(|Rvbicki fc Lightman|[l98^ is 



1 a/ 

c 'dt 



-npl + £ , 



(1) 



where I is the specific radiation intensity, n is the opac- 
ity coefficient, p is gas density and e is emissivity of the 
gas. In general I and e are functions of direction, radiation 
frequency, position and time. The first term on the right 



hand side represents removal of radiation from the beam by 
absorption and scattering, whereas the last term describes 
local emission of radiation into the beam's direction. 

In Monte Carlo methods, the radiation field is sam- 
pled via photon packets. In the simplest reincarnation of 
the method, both terms on the right hand side of equation 
[T]are treated stochastically. The photon mean free path, A, 
is calculated as A = l/{Kp). A random number, ^, uniformly 
distributed over [0, 1], is generated. A photon is allowed to 
travel distance Al — — Aln^, at which point it interacts 
with gas by passing its energy and momentum to gas com- 
pletely. The photon is then re-emitted according to physics 
of a chosen set of radiation processes, and followed again 
until it escapes the system. 

Our radiation transfer scheme is slightly modified from 
this. Firstly, instead of using discontinuous photon jumps, 
we explicitly track packet's trajectory in space as 



r(t) = ro 



(2) 



where r(t) and ro are the current and initial photon lo- 
cations, and |v-y| — v.y = const is the photon propagation 
speed (not necessarily the speed of light; see ^2.5|) . Secondly, 
photon momentum, p^, and photon energy E.y = cp^ are re- 
duced continuously due to absorption or scattering as 



dpj 



' A 



(3) 



Thus, the first term on the right hand side of equation 
[T] is treated continuously rather than stochastically. This 
reduces the statistical noise significantly. Photons are dis- 
carded when their momentum drops below a small (~ 10"'*) 
fraction of their "birth momentum" p^^o- 

However, the re-emission term, i.e., the last one in equa- 
tion [T] is modelled similarly to classical Monte-Carlo meth- 
ods in that it is stochastic in nature. This is unavoidable 
in Monte Carlo methods due to t he need to e mploy a finite 
number of photon packets (e.g., ILucvI [l999l ). In practice, 
at every photon's time step At-y, the photon momentum ab- 
sorbed, Ap-y, is calculated according to equation[3] As stated 
in the Introduction, we concentrate in this paper on the ra- 
diation pressure effects, and assume that radiation absorbed 
from the beam is completely re-emitted on the spot in a new 
random direction. This is the case for pure scattering of ra- 
diation or for local radiative equilibrium between radiation 
and gas. In these cases the amount of momentum (energy) to 
be re-emitted in the interaction considered is exactly Ap-,. 

Probability of re-emitting a new photon with momen- 
tum (energy) p^o (cp^o)is defined as — Ap-y/p^o- A ran- 
dom number, ^, uniform on [0, 1], is generated. If ^ < Wj, a 
new photon with momentum p^o and a random direction of 
propagation is created. This conserves photon field's energy 
and momentum in a time-average sense. 

This scheme can be adapted to allow for non- 
equilibrium situations when radiative cooling is not equal 
radiative heating. Multi-frequency radiation transfer is also 
straight forward if tedious to implement. We have done some 
of these developments already and will report it in a future 
paper. 
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Figure 1. Comparison of radiation energy density inside a slab of matter computed via our radiation transfer formalism (red dashed 
curves) and by the traditional Monte Carlo technique (solid). The difference is minor, and is entirely due to a finite time- resolution of 
the red curve, as explained in the text. 



2.2 A static test of radiation transfer 

The radiation transfer approach described in i]2.1l embodies 
the classical radiation transfer equations and must be correct 
on that basis. Nevertheless, we tested our approach against 
the discrete Monte- Carlo method on a simple test problem 
where the gas density and opacity are fixed and independent 
of the radiation field. This setup does not require using gas 
dynamics at all, and hence no SPH particles were used to 
model the density field. In the test, photons are emitted 
at a constant rate at the midplane of a plane-parallel slab 
z = [—1,1]. The main parameter of the problem is the total 
optical depth of the gas between z = and z = 1, Tt. We ran 
the radiation transfer part of the code until a steady state 
photon packet distribution was achieved. 

Figure [T] shows the ratio of the local radiation energy 
density to the radiation fiux out of the slab as a function of 
z-coordinate inside the slab. For convenience, the speed of 
light is set to unity. Two tests are shown, one for Tt = 2 (left 
panel) and the other for Tt — 20 (right panel). Black solid 
curves show the same quantity computed with the standard 
Monte Carlo technique, whereas the red dashed curves are 
computed with the modified method described above. For 
the latter time-dependent calculation, we averaged over a 
number of snapshots. The disagreement is quite small ev- 
erywhere except for the region near z = 0. The difference in 
that region is simply due to the finite time resolution of the 
red curves. Namely, these curves are obtained by averaging 
snapshots of the simulation data. The time step between 
the snapshots is long enough for the photons introduced at 
times between the snapshots to diffuse away from their start- 
ing z = location, broadening the peak. The Monte Carlo 
curve, on the other hand, is obtained by following individual 
photons' progress out of the slab, and therefore the informa- 
tion about the initial location of the photon is retained and 
is evident in the peak ai z — 0. These tests demonstrate 



that the radiation transfer is modelled properly with our 
approach at least in the simple situations considered here. 

2.3 Interactions between photon packets and SPH 
particles 

We now discuss coupling the radiation transfer scheme to 
SPH. A photon mean free path is given by 



where k is the gas opacity coefficient and p is gas mass den- 
sity at the photon's current location. We require an accurate 
method of local gas density estimation, and therefore we use 
directly the density field representation of the SPH method. 
In the SPH, gas density at a position r is given by the su- 
perposition of smoothed contributions from individual SPH 
particles: 

p{v) = Y.m,W{\r-n\,h,)^Y.p,{r) , (5) 

i i 

where W is the SPH kernel, ri is the location, hi is the 
smoothing radius of particle i, and pi = rmWi is the contri- 
bution of particle i to the density at r. The summation goes 
over all the "neighbours" - SPH particles that contribute 
to p at r, i.e., have non-zero values of W at this location. 
Note that in this approach there is no expectation of the 
number of gas neighbours of a photon, N^n, to be constant 
or limited in any other way. In particular, having A^gn = is 
perfectly feasible and simply means that the photon propa- 
gates through an "empty" patch of space. 

Having determined A at the photon's location, we find 
the decrement in photon's momentum, Ap-,, at that loca- 
tion according to equation [S] This decrement is then passed 
to the photon's SPH neighbours to enforce conservation of 
momentum. If Agn > 1, there is a question of what fraction 
of Ap-y should be passed to a particular neighbour i. We as- 
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sume that the neighbours contribute to the interaction with 
a photon directly proportionally to their density p; at the 
photon's location (see equation [5]) . Hence the momentum 
passed to neighbour i is 

Ap^. = n^Ap^ . (6) 
P(r) 

As we limit the photon time step to ensure that At = 
w-yAt-y/A <^ 1 (see equation [Qjl , this approach yields physi- 
cally correct result. 

The momentum transferred to the SPH particle i from 
interactions with different photon packets 7 is additive, and 
is used to define the radiation pressure force on that particle. 
We sum all the interactions that the SPH particle experi- 
enced during its time step At, and then define the radiation 
pressure force on particle i by 



frad,i 



T,y AP7« 



At 



(7) 



Due to limitations on photon propagation time steps 
At J (see below), an SPH particle may interact several times 
with a particular photon 7. There is however ho paradox 
here, eis this simply means that the radiation transfer equa- 
tion is integrated in multiple points along the photon's tra- 
jectory within each SPH particle, increasing the precision of 
the method over traditional Monte-Carlo approaches. This 
formulation also allows us to enforce an exact pair-wise con- 
servation of momentum in the interactions of matter and 
radiation. The energy transfer can be calculated in exactly 
same manner except Ap.^ is replaced by A^E.^. 



2.4 Photon packet creation and propagation 

External radiation sources such as stars or accreting com- 
pact objects emit photon packets with a given momentum 
p-yo- The rate of packet emission is given by 



Ny = 



cp-,0 



(8) 



where L is the luminosity of the source and c is the speed 
of light. The radiation field can be chosen to be isotropic or 
beamed/restricted to a range of directions. 

There are several issues to consider when deciding how 
far the photon can travel in a single fiight Al — v-yAt-f. First 
of all, this distance should be much smaller than the typical 
SPH smoothing length, h, in the region, or else the photon 
will "skip" interactions with some SPH particles altogether. 
This constraint is important in both optically thin and thick 
regimes. In the optically thick case, an additional constraint 
need to be placed to insure that photons do not propagate 
in one step by more than a fraction of their mean free path, 
A. In the opposite case photons would "diffuse" though op- 
tically thick regions in an unphysical way, i.e., too quickly. 
These constrains are combined by requiring 



Aty = iStmin 



A 



(9) 



where 5t <^ 1 is a small dimensionless number. In practice 
we use 5t = 0.03 — 0.3. 

As discussed below in H2.5I in the "prompt escape" 
regime, it is possible to reduce the photon propagation speed 
below the speed of light, which then allows us to integrate 



photon trajectory on longer time steps without compromis- 
ing the physics of the problem. 

2.5 Regimes of radiation transfer 

Following iKrumholz et al.1 (|2007l '). radiation hydrodynamics 
of a problem can be divided into three different regimes. 
Let u be a characteristic gas velocity, such as the sound 
speed or the bulk gas velocity, whichever is greater. Define 
/3 — u/c and optical depth of the system r — l/X, where / is 
the geometric size of th e system. The three limiting regimes 
IKrumholz et al.ll2007f ) are 

T <^ 1 the free streaming limit, (10) 
r 2> 1, /3r <C 1 the static diffusion limit, (11) 



r 2> 1, /3r 2> 1 the dynamic diffusion limit. 



(12) 



In the first regime a typical photon leaves the system 
in a single flight. In the second case the photon scatters or 
get absorbed and re-emitted approximately S> 1 times 
before leaving the system. For what follows it is important 
that in the first two regimes photons escape from the system 
on timescale, tcsc, much shorter than the matter distribution 
can alter significantly, i.e.. 



= - 1 +r < - 
c u 



(13) 



Here -R is the geometric size of the system. In contrast, in the 
dynamic diffusion limit the photon diffusion time is larger 
than the crossing time of the system, R/u. 

Equation [13] shows that the exact value of the speed 
of light is irrelevant in the free streaming and the static 
diffusion limits; it is so high it can be considered infinite. To 
the contrary, in the dynamic diffusion limit the exact value 
of speed of light is important as it defines the time scale on 
which radiation from the system leaks out, and that time 
scale is long. 

Therefore, we combin e the free streaming l imit and 
the static diffusion limit of IKrumholz et al] (|2007l ) into the 
"prompt escape" regime given by the equation 1131 In this 
limit it is numerically convenient and physically permissi- 
ble to reduce the photon propagation speed u-y below c, as 
long as the system still satisfies equation [13] with c replaced 
by v.y. We found that this speeds up calculations in which 
gas dynamics is important, although the scaling is not as 
efficient as the factor c/v.y. The reason for that is that al- 
though photon time step (equation [5| is indeed longer by 
the factor of c/u-y, the number of photon packets for a given 
p-y is correspondingly higher (see equation [8] and note that 

N-, ~ Nytesc). 



2.6 Implementation in Gadget 

This radiation-gas momentum transfer meth od has been im - 
plemented in the SPH/N-body code Gadget (|Springelll2005h 
that is widely used for cosmological simulations. For the 
tests presented here the cosmological options of the code are 
turned off. Gadget uses a Barnes-Hut tree to speed up cal- 
culation of gravitational forces and for finding neighbours. 
As photon packets have no mass associated with them, we 
turn off the gravity calculation for these particles. They are 
also not included in building the Barnes-Hut tree. 
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To find SPH neighbours of a photon, we first find aU 
SPH particles that are inside a sphere with size ftaearch which 
is chosen to be much larger than the mean SPH particle 
smoothing length. We then further select only those SPH 
particles i that contain the photon within their smoothing 
length hi. 

After calculating the radiation pressure force for a given 
SPH particle i, the corresponding radiative acceleration 
Srad.i = frad.i/'Tii is added to the hydrodynamical and gravi- 
tational accelerations that the particle experiences. We have 
also added the radiation pressure acceleration to the time 
step c riteria for the SPH particles, as described in ISpringell 
(|2005l ). This ensures that SPH particle time steps are ap- 
propriately short for particles with large radiation pressure 
accelerations. 




2.7 Static SPH radiation transfer test 

In !^2.2l we tested the radiation transfer methods for the den- 
sity field given by a simple analytical function (a constant). 
In i|2.3l we presented a way to model radiation transfer in 
an arbitrary density field represented by the SPH particles. 
It is a logical step forward in complexity of the tests to now 
repeat the slab test of i|2.2l in a self-consistent SPH density 
field. 

To accomplish this, we consider a non self-gravitating 
accretion disc in orbit around a massive (M = 1 in code 
units) central source. The disc is assumed to be locally 
isothermal, with internal energy u scaling as it = uo{Ro/R), 
where R is radius, and uo — 0.02 is the internal energy at 
the inner edge of the disc, Rq = 1. This scaling yields a 
constant ratio of the vertical disc scale height, H, to ra- 
dius. The outer radius of the disc is -Rout ~ 3. The disc 
was then relaxed for a large number of orbits without any 
radi ation field, which yields a Gaussian vertical density pro- 
file (jShakura fc Sunvaev|[l973 ) . Keeping this density profile 
fixed, we then introduced photon sources in the disc mid- 
plane and allowed the photons to propagate out of the disc 
keeping the opacity coefficient k fixed everywhere in the disc. 

Figure [2] shows the resulting photon energy density dis- 
tribution within the slab (solid curve) in exactly same man- 
ner as in the fixed density tests shown before in Figure [1] 
As before, the radiation field has been averaged over several 
snapshots to reduce statistical noise. The respective Monte- 
Carlo result obtained for the same density distribution is 
shown in Figure[2]with the red dashed curve. There is a very 
good agreement between the curves, except for the peak re- 
gion in the red curve. Time sampling differences in the two 
simulations explain the discrepancy in the curves; see discus- 
sion of Figure [T] in H2.2I This test suggests that the radiation 
transfer part of our approach is working as expected. Below 
we shall move on to tests that involve gas dynamics. 



3 ABSORPTION ON THE SPOT TESTS 

Young massive stars produce most of their radiation in the 
UV wavelengths. Dust in the interstellar gas has a very 
large absorption opacity for UV photons, with fcuv up to 
a few hundred cm^/ g. If these photons are absorbed and 
re-emitted in the infrared, where Roseland opacity is orders 
of magnitude smaller than kuv, in many applications one 
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Figure 2. The ratio of the radiation energy density to tlie radi- 
ation flux for the simulation described in section 12.71 (black solid 
curve) compared with a traditional Monte-Carlo calculation for a 
slab with same optical thickness (red dashed curve). While sim- 
ilar in spirit to tests shown in Figure [T] the present simulation 
uses an actual SPH density from a live accretion disc simulation. 

can effectively assume fcir = 0. In this case it is sufficient to 
consider only the radiation pressure from the UV photons 
emitted by the young stars, and neglect the re-radiated com- 
ponent. Finally, one can approximate the large UV opacity 
by an infinitely large one, fcuv = oo. Photon packet prop- 
agation is then trivial: packets travel in straight lines with 
constant momentum p-y until they encounter an SPH parti- 
cle(s), at which point they are absorbed and their momen- 
tum is transferred to that particle(s). 

This "on the spot" absorption method can also be used 
to model fast gas winds from massive stars or luminous black 
holes in the momentum driven regime. In the latter case 
the cooling time in the outfiow is short. Shocked outfiow 
gas cools very quickly, and hence its thermal pressure can 
be neglected. The momentum transferred to the ambient 
medium provides the push to drive the shell out. If the wind 
velocity is much higher than the velocity of the expanding 
shell and the speed of sound in the ambient gas, one can also 
neglect the mass outfiow from the source. This is justified 
as the mass fiux in the wind is small compared with that of 
the ambient gas being driven out. 

The "on the spot" approximation presents a convenient 
test ground of our code since it is possible to derive exact 
analytical solutions in the simplest cases. In particular, we 
consider a single radiation source embedded in an infinite 
initially uniform isothermal medium. Gravity is turned off 
for simplicity. 

In practice we set up periodic boundary conditions for 
a cubic box with dimensions, I, of unity on a side. These 
boundary conditions are appropriate and do not affect our 
results since the radiation force effects are contained to a 
small region within the box during the simulations. Gas in- 
ternal energy is fixed at u = 1, hence sound speed Cs = 1. 
The total mass of the gas inside the box is M = 1. The unit 
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of time is l/cs- 10® SPH particles is used in these tests. The 
initial condition is obtained by relaxing the box without the 
radiation source to a state of nearly constant gas density 
over hundreds of dynamical times for the simulation box. 



3.1 A steady-state case 

We first consider a case with a relatively small source lu- 
minosity L-y = 0.025 in the code units. In this case a quasi- 
steady state should be set up quickly. We consider that state 
here. Figure |3] shows the density profiles obtained in four 
different runs at dimensionless time t — 3.5, when the SPH 
quantities approach a quasi steady state. Theoretical ex- 
pected density profile is given by p = for 7? < 7?cav and 
p = const 1 for _R ^ Rcav Here -Rcav is the size of the 
cavity opened up by the radiation. This is obtained by re- 
quiring a force balance between the momentum flux from 
the source, L-y/c, and the external pressure of the gas pcxt- 



Rc 



1 1/2 



(14) 



The expected discontinuous density profile is shown with the 
solid (black) line (we have neglected a slight density increase 
due to evacuation of the gas from the cavity, since -Rcav ^ I, 
the box size). Figure|4]shows instantaneous velocities of SPH 
particles corresponding to the respective panels in Figure |3l 
Theoretically expected result for velocity in the steady state 
is u = everywhere, of course. 

Figures [3] and [3] show that at relatively large values 
of the photon packet momentum, corresponding to a lower 
packet injection rate A'^ — L^/{cp^), the method is inaccu- 
rate. For p-y = 5 X 10"'', no cavity is opened at all around 
the source. Moreover, SPH particle velocities have random 
components at a 30% fraction of the sound speed, even at 
regions far behind the expected discontinuity at i? = -Rcav 
In contrast, tests with p^i = 10~^ and 10~* produce nearly 
identical results with velocity fluctuations below a few per- 
cent level at 2-Rcav . 

Physically, if > Psph, interaction of an SPH particle 
with a single photon packet may accelerate the particle to 
velocity exceeding the sound speed, leading to a significant 
numerical noise. In the tests presented here, a typical SPH 
particle momentum is Psph = m-sphCs — 10~®, as msph = 
10~®. Therefore, the two runs with > psph do very poorly, 
as expected. These numerical experiments show that the 
minimum accuracy requirement in optically thick limit is 



P^ rvj Psph 



(15) 



For further quantitative analysis of the errors, consider 
now the width A-R of the transition layer. This layer is de- 
fined as the region in which SPH density goes from zero to 
unity. First, note that on average, a gas parcel within vol- 
ume ~ situated on the inner face of the transition region 
receives one photon packet every 



1 47r-R2^, 

iitl ~ — To — 

7V^ vr/i2 



(16) 



seconds. During this time, SPH particles inside the parcel 
experience no radiation pressure "pushes". We assume that 
in the absence of radiation forces, the volume is quickly ac- 
celerated by the pressure gradient force to i;_r ~ — Cs. There 



are ~ Nnb SPH particles within the volume, where Nnh is 
the typical number of SPH neighbours (usually chosen to be 
around 40). Due to the spherical symmetry of the problem, 
the parcel travels inwards to the radiation source a distance 
A-R = vitAti before it is turned back by the arrival of the 
next photon packet. Finally, since the volume element is lo- 
cated on the inner boundary of the cavity evacuated by the 
radiation pressure, it will expand when moving into the cav- 
ity. Hence the appropriate smoothing length of the element 
will be larger than that far from the cavity, and in general 
we expect h ~ AT? in that region. Using the latter estimate 
in equation 1161 we obtain for the thickness of the transition 
layer. 



AR 

2-/?cav 



1/3 



(17) 



The transition layer thickness calculated in this way is in- 
dicated in Figure [3] with dotted (-Rcav — A-Rtr) and dashed 
(-Rcav -I- A-Rtr) lines. It appears to be a good estimate for the 
three runs out of the four shown in the Figure. In the run 
with the highest photon packet resolution, i.e., the smallest 
Pt, equation [17] is an under-estimate. This is simply because 
a density discontinuity in SPH cannot be narrower than the 
width of the kernel, i.e., the minimum smoothing length, 
which is about h — 0.02 for all the four tests. Inserting equa- 
tions [14] and [8] into equation 1171 and requiring A-R to not 
exceed the characteristic smoothing length of the problem 
leads to the same requirement of the packet's momentum to 
be small enough compared with the typical SPH momentum 
(equation [15]). 

Summarising the results of these tests, we see that we 
can use the condition [15] as the photon packet resolution 
requirement in an optically thick limit. In addition, stochas- 
tic noise arguments appeared useful when determining the 
width of the transition layer in Figure [3] 



3.2 Momentum-driven wind test 

Keeping to the model of a single source radiating photons 
isotropically into an infinite, constant density, isothermal 
medium, we performed a further test whereby the luminosity 
of the source was significantly higher, = 50. This time 
the focus of our attention is the evolution of the radius of 
the expanding shell with time, well before it reaches the 
equilibrium cavity size. 

This test is closely related to the ram-pressure driven 
wind test (e.g., IVishniad ll983l : iGarcia-Segura et al ] ll996l : 
[Dale fc Bonnelll l2008l ). but with the important distinction 
that in our method the photon 'gas' is massless. This how- 
ever should be a good approximation to a low density but 
high velocity outflow. 

The other important aspect of our model is the isother- 
mal condition, corresponding physically to a shocked gas 
region that is allowed to cool very quickly, i.e. when the 
cooling time in the shocked gas is far shorter than the dy- 
namical time. This tends to occur when the temperature 
is less than lO'' K, so that line cooling from ions domi- 
nates (jLamers fc Cassinellilll999l ). The picture then is that 
of a momentum-conserving 'snowplow' phase of an expand- 
ing gaseous bubble as it sweeps up the ambient gas into an 
increasingly massive shell. This model has three zones: (i) 
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Figure 3. Density of selected SPH particles versus radial distance from the source for tests described in f[3] The dimensionless luminosity 
of the source is L-y = 0.025 for all the tests. The momentum of the individual photon packets varies between p-, = 5 X 10~^ to p-y = 10~*, 
as indicated in each of the panels. The solid line showrs theoretically expected density profile, whereas the dotted and the dashed lines 
show the estimated error. 
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Figure 4. Absolute velocity of SPH particles for tests shown in Figure [3] Theoretically expected solution would have i; = everywhere. 
As before, the results in the two upper panels, corresponding to the cases of a largo photon momentum, arc quite inaccurate and show 
significant numerical fluctuations in v. The lower panels show such fluctuations only in the narrow layer adjacent to the inner cavity. 
Fluctuations behind the layer are very subsonic and can be neglected for practical purposes. 
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an evacuated cavity through which free-streaming photons 
are being emitted from the central star (ii) a shell of shocked 
ISM gas at temperature To bounded by an outer shock front 
and (iii) the undisturbed uniform ISM at temperature Tq. 

The initial conditions for this test and numerical setup 
are same as in i)3.1l The choice of the higher luminosity is 
motivated by a requirement that the Mach number be high 
so that the pressure from the ambient gas did not have a 
significant efi'ect on the evolution of the shell during the 
run. To aid this, the gas internal energy was set to a lower 
value of M = 0.1. 

The analytical solution to a momentum-conserving bub- 
ble can be derived by considering the thin-shell approxima- 
tion, whereby the swept-up gas is assumed to be concen- 
trated in an infinitely thin shell that is being driven by the 
impinging wind. In reality the shell gets thicker the more 
gas it sweeps up, but at all time s the thickness of the she ll, 
D, is much less than the radius jClarke fc Carswell|[2003l ): 



D 



R 
3M2 



(18) 



where S> 1 is the Mach number and R is the radius of 
the shell. 

The zeroth order accuracy analytical solution is found 
by equating the rate of change of momentum of the expand- 
ing shell to the momentum flux from the stellar luminosity: 



dt 



(^^■kR^Po]R 



L 

c 



(19) 



The time evolution of the radius of the shell is therefore 
given by. 



R 



( 3L 



y l-KCpQ 



1/4 



,1/2 



(20) 



We assumed here that the velocity of the particles carrying 
the momentum outflow from the source, ii-y, is much larger 
than the shell velocity, 7?. 

This solution can be further improved by taking account of 
the finite speed of the wind (photon packets in our code). In 
this case there is a time delay of Rjv^ before the wind par- 
ticles strike the ambient gas. In addition, external pressure 
from the ambient medium provides some restoring force, 
slowing down the expansion of the bubble. Factoring these 
effects into equation [T5] leads to a more accurate description 
of the shell: 



R 



R 



4nR^ pc^t 



(21) 



which is integrated numerically to obtain a solution for R. 

Figure [5] shows the results of the test with p-y — 10^^. 
The solid black line shows the time evolution of the shock 
front (identified by the peak in the density distribution) in 
the SPH simulation. The analytical result given by equation 
I20l is shown with the red dotted line. The green dotted line in 
Figure [5] is a numerical solution of the more accurate equa- 
tion [21] The difference between the two expected solutions 
(red and green) is rather minimal since the external pressure 
is small compared with the ram pressure of the outfiow and 



the photon packet's velocity is large, v-y 



200. 



Figure [5] shows an excellent agreement between the ex- 
pected ID solution and the simulation's result. 




0.005 



0.025 



Figure 5. Radius of a shell expanding into a uniform medium as 
a function of time (see ^3.21 The shell is driven by an isotropic 
outflow from a point source. The solid black line is the evolution 
of the density peak (corresponding to the position of the shock 
front), the red dotted line is the analytical solution without cor- 
rections as per equation [19] and the dotted green line is the more 
accurate analytical solution corrected for a finite photon speed 
< c and the pressure of the external medium. 



4 PLANE-PARALLEL SLAB TESTS 

A simple but useful test of gas dynamics is provided by 
optically thin gas illuminated by a constant radiation fiux. 
Clearly the radiation force in this case must be a constant 
independent of position or time. In the stochastic approach 
of our method, however, the force can in principle vary spu- 
riously due to statistical fiuctuations, which may lead to nu- 
merical artefacts. Our goal here is to investigate their nature, 
magnitude and scaling with the number of photon packets 
used. 

As in ^ the initial condition is a uniform density 
box with a side equal to 1. The SPH particles are al- 
lowed to interact with each other hydrodynamically only 
(i.e., self-gravity is turned off). In the tests reported here, a 
given number of photons A^.^ propagates in the positive x- 
direction. The momentum of photons is set by requiring the 
time and volume averaged radiative acceleration be equal to 
a parameter ao- 

In the limit of an infinite number of photons and an 
absence of any numerical artefacts, all the SPH particles 
should be accelerated with acceleration ao, so that their x- 
component of velocity is Vx = aot. In reality different SPH 
particles interact with different number of photons due to 
stochastic fiuctuations, and hence their accelerations and 
Vx are not exactly equal. We define the Mach number of 
spurious SPH velocity fiuctuations as a measure of the error 
in this test: 



Morr(t) = — 
Cs 



-5:(v>-v(t))^ 



1/2 



(22) 



where A'' is the total number of the SPH particles, and v{t) = 
Vx{t) = aot is the average particle velocity at time t. The 
quantity in the square brackets is the velocity dispersion, of 
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We can build a simple "theory" to estimate the er- 
ror in these tests. We shall assume that velocity dispersion 
is driven by the stochastic fluctuations in the number of 
photon-SPH interactions in different regions of the box. We 
shall also assume that fluctuations in these will be uncor- 
related on a timescale of order the sound crossing time of 
the box, Lbox/cs- During this time the mean number of the 
interactions between photons and a parcel of gas with size 
(where h is the smoothing length) is 



^box 



(23) 



where — N-y/L^^^ is the volume average density of pho- 
tons in the box. The parcel is going to be accelerated to 
velocity aoLbox/cs during this time. Each photon then ac- 
celerates the parcel by a velocity increment of 



Ai;i = ao-Lbox/cs/A''pa 



(24) 



As the number of photon-parcel interactions fluctuates by 



= iV, 



1/2 



we expect that random velocity fluctua- 



tions will be of the order of 



Aw = /\viNI 



1/2 

ass 



(20 -Lb 



1/2 



(25) 



Figure |6] compares the measured error and the one predicted 
by equation [25] for three runs. In the runs the SPH particle 
number is fixed at A'sph = 10'', ao — 0.1, = 10, and the 
number of photon packets inside the box is varied from 30 
(upper curve), through 300 (middle, dashed curve) to A^'-y = 
3 X 10"^ (lower, dashed curve). The respective straight lines 
show the "theoretical" estimate given by equation 1251 The 
agreement is reassuringly good. One can also note that the 
error indeed saturates rather than grows with time, except 
for a few initial sound crossing times. This is due to the 
stochastic fluctuation becoming uncorrelated on dynamical 
time, as we assumed above. 

We also ran a large number of additional tests all with 
the same setup described in this section, but varying the 
number of SPH particles up to 2 x 10^, mean radiation pres- 
sure acceleration ao (increasing it to 1 and 10), and also 
varying the photon velocity v-y by a factor of 2 in either 
direction. The resulting errors and their scalings were com- 
pared with that predicted by equation 1251 conflrming the 
excellent agreement of the expected and measured errors 
further. 

Equation 1251 shows that, at the same ao (equivalently, 
the same radiation flux), the velocity error scales as 



Ati cx 



1/3 



sph 



hN. 



1/2 



A. 



1/2 



(26) 



The scaling of Aw with the number of photon packets as oc 

— 1/2 

A-y could be intuitively expected on the basis of Poisson 
statistics of random fluctuations. The scaling of equation 
1251 with the smoothing length h, which is proportional to 
Agpj^ , shows another important point. If a larger number of 
SPH particles is used, resulting in a higher spatial resolution 
within the simulation volume, then the number of photon 
packets should also be increased (if one wishes to keep the 
velocity errors within a given limit). 
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Figure 6. Velocity error defined as the Mach number of velocity 
fluctuations (equation [22]l versus time for three optically thin slab 
tests described in ^ The number of photon packets is varied 
from 30 (upper thick solid curve), through 300 (middle, dotted 
curve) to A'-y = 3 X 10"^ (lower, dashed curve). The respective 
straight lines show the "theoretical" estimate of the error given 
by equation 1251 



5 A SUSPENDED OPTICALLY THIN CLOUD 

Consider optically thin gas in the vicinity of a very massive 
point mass, e.g., a black hole, radiating exactly at the Ed- 
dington limit. Obviously, gas particles should experience no 
net force from the point source. In particular, if SPH parti- 
cle velocities are zero initially, they should remain zero, and 
any gas motion is an evidence of problems in the numerical 
methods used. 

For ease of analysis, we consider a self-gravitating gas 
cloud, with a polytrophic index F = 5/3. In the absence 
of external forces the cloud settles into a well known equi- 
librium configuration for polytropic stars. This equilibrium 
state is obtained by relaxing (that is evolving) the isolated 
cloud for many dynamical times to provide a noise-free ini- 
tial condition. 

Gravity is switched on in this test. However, as our focus 
is on testing the radiation pressure force, the gravitational 
force of the gas acting on the point mass is switched off. The 
radiation source is hence fixed in its initial location, which 
simplifies a quantitative analysis of the results below. 

The mass unit for the test is the mass of the point 
source, M. The mass of the cloud is 0.023. The polytropic 
constant is such that the equilibrium size of the cloud, Rc\, 
is approximately 4 units of length. 25,000 SPH particles are 
used to model the cloud. 

The tidal force from the central point is of the order of 
Ft ~ 2(GMMc\/B?)Ra, where R is distance to the point 
mass. There are two distinct regimes, then. In the first the 
tidal force acting on the cloud exceeds the self-gravity of the 
cloud, _Fsg ~ [GM^i/ R^i), whereas in the opposite regime the 
self-gravity of the cloud dominates over the tidal force. The 
former regime occurs when R — zo, the separation between 
the point mass and the centre of the cloud, is relatively small. 
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zo ^ Rci- The case of a negligible tidal force, Ft <§; Fsg, 
occurs when the cloud is far away from the black hole, e.g., 

5.1 z = 10 tests 

In these tests the cloud is initially positioned at 2; = and 
the emitter is at 2; = —10. At this separation, the tidal force 
is approximately 10 times greater than self-gravity of the 
cloud. Thus, a ~ 10% error in the radiation pressure force is 
larger than the self- gravity holding the cloud together. These 
tests are hence expected to be very sensitive to numerical 
deficiencies of our method. 

We performed two tests in this setup. In the first, la- 
belled LRZIO (low resolution, z = 10), the dimensionless 
photon momentum is p.^ = 2 x 10~*, whereas in the sec- 
ond, named HRZIO (higher resolution), p-y = 2 x 10~^. 
The mass of the SPH particles is the same in both tests, 
mspH f« 10~^. The typical momentum of SPH particles, de- 
fined as psPH = mspH^yGM/zo, is psph « 4 x 10~^ in code 
units. Therefore, in both of these runs p-y <ti Psph. 

Figure [7] shows the column density profile of the cloud 
at time t = 100 for both runs. The left panel shows run 
LRZIO and the right panel shows HRZIO. Time t = 100 cor- 
responds to about 3 free-fall times at the given source-cloud 
separation. In both cases there is a certain deformation of 
the cloud. As expected, the lower resolution test produces 
poorer results than the higher resolution test. A more de- 
tailed error analysis is described in i)5. 31 below. 

5.2 z=30 test 

In this test the separation of the cloud and the source is set 
to Zo — 30, thus we label the test Z30. The photon packet 
momentum is pj = 10~* in the test. The parameters of the 
gas cloud are exactly the same as in tests LRZIO and HRZIO, 
but because the cloud is farther away from the point mass, 
typical SPH momentum is a little lower, pspH ~ 2.5 x 10~^. 
Thus in terms of ratio p^/psPH test Z30 is very similar to 
LRZIO. 

However, since the cloud is farther away, the self-gravity 
of the cloud actually exceeds the tidal force by a factor of 
about 3. Therefore, one expects intuitively that the same 
level of stochastic error in the radiation pressure force will re- 
sult in smaller perturbations as the radiation pressure force 
itself is smaller. This is indeed borne out by the test. By 
the time t = 500, which is again about three free-fall time 
scales from the clouds' initial location, no visible deforma- 
tion of the cloud has occurred. We therefore do not show 
the column density for test Z30. 

5.3 Error analysis 

To quantify the accuracy of the tests, we calculate the centre 
of mass velocity of the cloud, function of time. 

Using this in place of v(t) in equation 1221 we define the 
mean Mach number of spurious velocity fluctuations in the 
cloud. Figure |8] shows these quantities as a function of time 
for LRZIO (red curves) and Z30 (black) runs. The Mach 
number of velocity fluctuations are shown with thick solid 
lines. The centre of mass velocities are shown normalised to 



0.3 F 




E 



c "■ : 

I -0.1 ^ ....... ..-.^ 

-0.2 F - 

100 200 300 400 500 
time 

Figure 8. The Mach number of numerical fluctuations (thick 
solid curves) and the average SPH velocity in units of 0.1 free- 
fall velocity (dotted) for tests HRZIO (red curves) and test Z30 
(black curves). Theoretically estimated Mach number errors are 
also shown with thin solid curves of same colours. 



0.1«K, where hk = \/ GM/zo, the circular Keplerian velocity 
at the appropriate separation zo for the test. 

The motion of the cloud's centre of mass characterises 
the systematic error of the method. Both runs demonstrate 
that the velocity of the centre of mass of the cloud is be- 
tween 1 to 2 % of the Keplerian velocity after about 3 free- 
fall times. This translates into inaccuracy in the calculation 
of the radiation pressure force versus point mass gravity, 
averaged over all SPH particles, of a fraction of a percent. 
This accuracy should be sufficient for simulations in which 
gravity forces are calculated to a similar fractional preci- 
sion. Furthermore, we found that this offset, although always 
small for a large enough number of photon packets, actually 
depends on the random number generator used to approxi- 
mate the isotropic photon field. The results can thus be im- 
proved by using quasi-random rather than pseudo-random 
sequences for photon packet's angular distribution. We shall 
investigate this issue in the future. 

The Mach number of the fiuctuations, instead, quanti- 
fies the difference in the radiative accelerations received by 
different parts of the cloud. We again use the simple logic 
of random stochastic ffuctuations in the number of photon 
packets A'paas interacting with an SPH particle with smooth- 
ing length h, as explained in ^ Thin solid curves in Figure 
[8] show the Mach number of velocity errors predicted by this 
simple argument. Evidently, the errors are again explainable 
by the Poisson noise estimates. 



6 A SPHERICALLY SYMMETRIC 
ACCRETION TEST 

In this test the initial condition is the same optically thin 
cloud as used in ^ but the point mass source is located 
exactly at the centre of the cloud, and allowed to accrete gas 
particles that are separated from it by less than R < Race = 
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Figure 7. Column density of gas in the tests LRZIO (left panel) and HRZIO (right panel). The optically thin cloud in the tests experiences 
gravity from the central source at z = — 10 and the radiation pressure force at the Eddington limit. The two forces should exactly balance 
each other. The snapshots are for dimensionless time t = 100, roughly three free-fall times. The higher resolution test HRZIO (lower p-y) 
shows a smaller degree of cloud deformation, as expected. 



0.2. The radiation field in this test is halved to L — 0.5 I/Edd. 
While the analytical solution to this problem (for the given 
initial density profile) is not kno-wn to us, we note that the 
setup is exactly identical to a non-radiating sink particle 
with mass A'l = 1/2 placed into the centre of the cloud 
and allowed to accrete the gas within the sink radius i?siiik- 
The desired control result is thus obtained by running the 
hydro-gravity part of the code only with the sink particle 
mass M = 1/2. Photon momentum is set to = 10~*. 

Figure [9] presents both the radiation hydrodynamical 
simulation and the control hydro run with M = 1/2 at di- 
mensionless time t — 3.5. The lower panel shows the SPH 
particle radial velocities normalised by the local free-fall ve- 
locity Vff = yjGM/R with M = 1/2. The upper panel shows 
the density profile of the cloud. Both velocity and density 
profiles of the RHD simulation and the control run are nearly 
identical, although one can notice a larger amount of scatter 
in the L = 0.5 I/Edd test. The results continue to be very 
similar at latter times as well, and the accretion rate his- 
tories are also very similar. The radiation transfer method 
thus performs reasonably well. 

To quantify the errors of the tests, we estimate their 
magnitude based on the Poisson noise arguments in the 
same manner as in ^ In particular, we can apply equa- 
tions [24] and 1251 except that the characteristic time is 
now not I/box /cs but the local dynamical time, t^yn{R) = 
R^'^/[GMY/^ with M = 1/2 (and G = 1 in the code 
units). Thus, the number of photon packets passed through 
a given particle is estimated as Afpass(-R) = Njtdyn{h/2R)'^ . 
The stochastic velocity fluctuations due to the finite number 
of photon packets is then 



Av ^ AviN^iliR) . (27) 

To compare the expected fluctuations given by equation 
[27] with those actually incurred, we ran a poorer resolution 
test with — 10~^, i.e., ten times larger than the test 
shown in Figure [9] The errors (scatter) are then obviously 
larger. Figure [10] shows the results in black dots for two 
snapshots at times indicated. The red lines show the control 
run, as before, but this time velocity curves in the lower 
panel of the figure include the expected errors. In particular, 
the upper red sequence of dots is Wcontroi + Av, whereas the 
lower red dots is iicontroi — Av, where Av is calculated with 
equation 1271 

Evidently, the scatter in the radiation transfer simula- 
tion (black dots) is very similar in magnitude to what is pre- 
dicted by equation 1271 At the earlier time snapshot, shown 
on the left side of Figure [Till the measured velocity scatter is 
somewhat larger than our prediction. We believe this is due 
to the fact that equation [27] does not include the Poisson 
noise resulting from the initial distribution of photons pack- 
ets at t = 0. The latter is generated by sampling a uniform 
random distribution in both emission time and directions. 

In conclusion, these tests demonstrate the potential of 
the method for accretion and radiation pressure problems, 
and that equation [27] once again provides a reliable estimate 
of the stochastic errors of the results. 



7 DISCUSSION 

A new method for implementation of radiation pressure 
force in an SPH code has been presented. Radiation is mod- 
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Figure 10. Similar to Figurc|9l but for a test with = 10 , i.e., lower resolution in photon packets. The red dots in the lower panels 
indicate the estimated range of errors (see text in §sec:accretion), which agrees reasonably well with the spread in the black dots. 



elled via a time-dependent Monte-Carlo approach. Photon 
packet's momentum is transferred to SPH particles via di- 
rect photon-to-gas interactions that are calculated as pho- 
tons propagate through the field of SPH particles. The SPH 
density field is used to calculate the SPH density at pho- 
ton packet's positions, thus maximising the accuracy of 
the method. As a result, local and global momentum con- 
servation in the interactions between radiation and gas is 
achieved. 

Several tests have been presented to check the consis- 
tency and accuracy of the method, as well as error scalings 
with parameters of the problem. It was found that the main 
deficiency of the method, as for any Monte Carlo method, is 
the Poisson noise resulting from a finite number of photon 
packets used. In order to reduce velocity fluctuations to sub- 
sonic levels, photon packet's momentum p-y should be chosen 
to be below the typical SPH particle momentum, mgphWsph, 
where Vsph is the characteristic SPH particle velocity. 

We shall now try to generalise and summarise the er- 
ror estimates and various constraints on performance of the 
code. 



7.1 Velocity errors 

We first assume an optically thin case. Let J-^ad be the radia- 
tion fiux at a given location in the gas. The fiux is related to 
the photon number density n^., photon momentum p-, and 
photon velocity v^y through 



rad ~ n^V^p^C , 



(28) 



A^pass ~ n^v^Txh 



nh' 



(29) 



The momentum passed to the SPH particle by one photon 
packet is 



''sph 



(30) 



and thus the velocity change Aiii = Kp^,/{Tvh^) is indepen- 
dent of SPH mass mgph (in the optically thin approxima- 
tion). Therefore, at time t the velocity fiuctuations can be 
estimated as 



Av{t) ~ Avi [iVpasst] 



1/2 



1/2 



(31) 



Now, let us define radiation acceleration time trad = 
't'sph/tSrad, where Orad = ^rad /"isph is the radiative acceler- 
ation acting on the SPH particle, and Frad = rn.sphnJ'ra.d/c 
is the radiation pressure force on the SPH particle. During 
this characteristic time, the radiation pressure force would 
accelerate the SPH particle from rest to velocity v ~ Vsph in 
the absence of other forces. 

Using eauation l31l the Mach number of the fiuctuations, 
defined by equation 1221 can be shown to be 

1 1/2 



Morr(i) 



Av{t) 

Vsph 



np^ 



t 



(32) 



Finally, introducing the "optical depth" of the SPH particle 
as Tsph ~ Kmsph/('7r/i^), we arrive at 

1/2 



A/orr(t) 



' s p h ~ 



Psph 



t 

^rad 



1/2 



(33) 



as the energy carried by the photon packet is p-yC. The rate 
at which photons are passing through SPH particles is 



where Paph = nispht'sph- This was derived in the optically 
thin limit, i.e., when Tsph <C 1. In the limit Tsph > 1, the 
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Figure 9. Spherically symmetric accretion test on the central 
source of mass M = 1 emitting at exactly 1/2 of the Eddington 
limit. The upper and the lower panels show the profiles of the 
SPH particle density and the ratio of radial velocity to the local 
free-fall velocity, respectively. The black dots show the results 
of the radiation transfer L = 1/2Le(J;j simulations, whereas the 
red ones show the results of accretion onto a non-radiating point 
mass with mass equal to 1/2. The two tests should yield identical 
results. 



momentum passed from one photon to an SPH particle is 
reduced from the expression given by equation 130 1 to 



(34) 



where A'nb is the number of SPH neighbours for the photon. 
As the latter is always at least one when a photon packet 
interacts with gas, we arrive at the following estimate for 
the optically thick case. 



Pi' 
Psph 



1/2 



t 

^rad 



1/2 



(35) 



Both the optically thin limit (equation [33} and the op- 
tically thick limit feauation l35p demonstrate the importance 
of choosing photon packet momentum, to be significantly 
smaller than that of an SPH particle to guarantee a reason- 
able accuracy. 

We also notice that applicability of the method to a 
particular problem depends on the desired level of accuracy. 
For example, tidal disruption of a gaseous cloud near a lumi- 
nous super-massive black hole is a highly dynamic process. 
To get an insight in the overall dynamics of the process, it is 
sufficient to simulate the system for a few dynamical times 
at the peri-centre of the orbit. A sufficiently small photon 
packet's momentum, p-y ~ 0.01 — O.lpsph, should provide a 



sufficient accuracy in this case. However if one is interested 
in a secular evolution of a system, then the required pre- 
cision is much higher, and careful tests should be done to 
establish the required value of p-y, or equivalently the photon 
injection rate A''.^. 



7.2 Number of photon packets 

Let us now estimate the required number of photon packets 
for a given luminosity of the system Lj, subject to the con- 
straint p^ — psph/r, r 3> 1. Number of packets emitted per 
unit time is 



CP7 



ph 



-A</gas ^sph 



(36) 



where we exploited the fact that Psph ~ (Mgas/Asph)fsph. 
Mgaa here is the total gas mass of the system. The minimum 
time that the simulation should last for is the dynamical 
time, R/vsph- During this time the number of photon packets 
emitted is N.y = N-yR/vsph, and thus 



= r 



L,{R/c) 



sph 



= r- 



sph "sph 



ph 



(37) 



where E^y = L-yR/c is the radiation energy emitted during 
the system light crossing time and Esph = MgasWgph is the 
total gas energy. If there are no additional sources of radia- 
tion in the cloud, such as bright stars, then the luminosity 
of the system is of the order of JSsph/icooi where tcooi 

is the cooling time. In that case the ratio is 



sph 



FR 

C^cool 



(38) 



Many astrophysically interesting problems are in the regime 
where the cooling time is comparable with dynamical time 
of the system, in which case the ratio becomes 



ph 



sph 



(39) 



From these expressions it is obvious that the number of pho- 
tons does not need to exceed the number of sph particles 
while satisfying the Merr 1 condition for these particular 
cases. However, if there is a very bright point source dom- 
inating the luminosity of the system, then the number of 
photon packets needs to be correspondingly higher. For ex- 
ample, in the tests presented in Sj5] photons produced during 
the duration of the simulations outnumbered SPH particles 
by factors of 10 - 100. 



7.3 Optical depth of the system 

Our radiation transfer method unfortunately becomes very 
expensive in optically thick medium, i.e., when the optical 
depth of the system, r, is very high. This is due to two fac- 
tors. First of all, photons scatter or get absorbed/re-emitted 
~ (that is many) times before they exit the system. They 
thus "hang around" for longer before exiting the system, 
hence increasing computational cost of the simulation. Sec- 
ondly, they require smaller time steps as we now show. 

Recall that in H2.5l we argued that it is numerically per- 
missible and beneficial to use photon packets with velocity 
smaller than the velocity of light in the following case. 
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In moderately optically thick, non relativistic plasmas, pho- 
tons pass through the system much quicker than dynamical 
time, and hence exact value of v-, is unimportant. The pho- 
ton time step can then be increased inversely proportional 
to w-y. However, this is possible only as long as 



V-y > Vsph (1 + t) , 



(40) 



or else using too small a value for Vj would incorrectly put 
the radiation transfer into the dynamic diffusion limit (see 
^2.5|l . Further, in the optically thick case the photon time 
step is limited to Atj = SfX/v^, where A is the mean free 
path (equation [9|. Comparing this to dynamical time of the 
system. 



At-, 

tdyn 



StX Vsph 

R v-, 



r(l + r) 



(41) 



where we used equation |40] to constrain v-,. This expression 
changes to 



At-y 



StV; 



'sph 



dyn 



(42) 



if one sets v-f = c. Obviously, in both of these cases the 
time step becomes very small as r increases, and the code 
becomes quite inefficient. 

In problems with a very high optical depth combining 
the dynamic Monte Carlo and the diffusion approximation 
in very optically thick regions would be optimal, although 
it is not clear whether such "on the fly" method could be 
devised. 



8 CONCLUSION 

A new time-dependent algorithm to model radiation trans- 
fer in SPH simulations was presented. In the present paper 
we concentrated on the radiation pressure effects, i.e., the 
momentum transfer between the radiation and the gas, as- 
suming a given equation of state for the gas. We performed 
a number of tests of the code. 

Our method is grid-less and can be applied to arbitrary 
geometries. The main disadvantage of the method, as with 
any photon packet based schemes, is the Poisson noise due 
to a finite number of photons. On the other hand, these 
stochastic errors can be readily estimated and controlled 
by increasing the number of photon packets. The method 
is best applicable to optically thin or moderately optically 
thick systems, as following photon trajectories become very 
expensive at high optical depths. 

Extension of the approach to include energy exchange 
between the gas and radiation will be presented in a future 
paper. Multi-frequency radiation transfer can also be natu- 
rally added. 
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